• We will use functions from packages base, utils, and stats (pre-installed and pre-loaded)
• We may also use the packages below (specifying package::function for clarity).

# Load pckgs for this R session

# --- General 
library(here)     # tools find your project's files, based on working directory
library(dplyr)    # A Grammar of Data Manipulation
library(skimr)    # Compact and Flexible Summaries of Data
library(magrittr) # A Forward-Pipe Operator for R 
library(readr)    # A Forward-Pipe Operator for R 
library(tidyr)    # Tidy Messy Data
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax

# ---Plotting & data visualization
library(ggplot2)      # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggfortify)     # Data Visualization Tools for Statistical Analysis Results
library(scatterplot3d) # 3D Scatter Plot

# --- Statistics
library(MASS)       # Support Functions and Datasets for Venables and Ripley's MASS
library(factoextra) # Extract and Visualize the Results of Multivariate Data Analyses
library(FactoMineR) # Multivariate Exploratory Data Analysis and Data Mining
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests
library(car)        # Companion to Applied Regression
library(ROCR)       # Visualizing the Performance of Scoring Classifiers

# --- Tidymodels (meta package)
library(rsample)    # General Resampling Infrastructure  
library(broom)      # Convert Statistical Objects into Tidy Tibbles

• Logistic regression is a type of Generalized Linear Model (GLM): a more flexible version of linear regression that can work also for categorical response variables or count data (e.g. poisson regression).

• Logistic regression is used to model a binary response variable: e.g. yes|no (or benign|malignant in the biopsy dataset).

• Logistic regression fits the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) to the data using the Maximum Likelihood Estimation method (unlike linear regression, which uses the Least Squares Estimation method).

If we have predictor variables like \(x_{1,i}\), \(x_{2,i}\), …, \(x_{k,i}\) and a binary response variable \(y_i\) (where \(y_i = 0\) or \(y_i = 1\)), we need a “special” function (or link function) able to transform the expected value of the response variable into the outcome we’re trying to predict.

If \(p_i\) is the probability that \(y_i = 1\)… \[ logit (p_i) = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \cdots + \beta_k x_{k,i} \]

…the logit function is defined as:

\[ logit(p_i) = \ln\left( \frac{p_i}{1-p_i} \right) \]

Name: heart_data.csv
Documentation: Toy dataset prepared for teaching purposes. See reference on the data here Data Analytics with R
Sampling details: This dataset contains 10,000 observations on 4 variables.

# Use `here` in specifying all the subfolders AFTER the working directory 
heart_data <- read.csv(file = here::here("practice", "data_input", "05_datasets",
                                      "heart_data.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

We could technically use linear regression (Figure 1) to study likelihood of having Heart Disease in relation to fast food spending, but it would not work:

• How should we interpret a continuous prediction when it is different from \(0\) or \(1\) ?
• Also, here the assumptions for linear model are not met:
  • Linear relationship between Y and X1
    
  • Residuals must be 1) approximately normally distributed and 2) uncorrelated
  • Homoscedasticity: residuals should have constant variance
  • Non collinearity: predictor variables should not be highly correlated with each other

# Hypothetical linear regression model
lin_reg_model <- lm(heart_disease ~ fast_food_spend, data = heart_data)
# Plot
ggplot(heart_data) +
  geom_point(aes(fast_food_spend, heart_disease, color = factor(heart_disease)), 
             alpha = 0.25, size = 1.5) + 
  scale_color_manual(values = c("0" = "#005ca1", "1" = "#9b2339")) +
  # regression line
  geom_abline(slope = coef(lin_reg_model)[2], 
              intercept = coef(lin_reg_model)[1], color = "black", size = 1.2, linetype = "dashed" ) + 
  labs(title = "Heart disease [0,1] against fast food spend (US$), with regression line") + 
  ylab("Heart Disease (Y/N)") + 
  xlab("Annual Fast Food Spend") 

• We can fit instead a logistic regression model to the heart_train data using:
  • the glm function for Generalized Linear Models;
    • with argument family = binomial to specify the logistic regression model which will use logit as the link function;
    • and after ~ we specify the predictor variables.

# Fit a logistic regression model
heart_model <- glm(heart_disease ~ coffee_drinker + fast_food_spend + income,
                   data = heart_train, 
                   family = binomial(link = "logit"))

• Table 1 (next) shows the model summary, with the coefficient estimate for each predictor.
• the broom::tidy function converts the model summary into a data frame.

heart_model_coef <- broom::tidy(heart_model) %>% 
  mutate('signif. lev.' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  # add odds ratio
  mutate(odds_ratio = exp(estimate)) %>%
  relocate(odds_ratio, .after = estimate) %>%
  # for visibility
  #select( -statistic, -p.value) %>%
  knitr::kable() %>% 
  #reduce font size
  kable_styling(font_size = 20)

heart_model_coef

• Intercept: the estimated intercept is typically not of interest
• Coffee Drinker: Coffee drinking has an estimated coefficient of -0.7296 which means that the coffee drinking is associated with a decrease in the odds that a person has heart disease. Specifically, it decreases the odds of heart disease (relative to non-coffee drinkers) by a multiplicative factor of e-0.7296, or 0.48, holding all other predictors fixed. Therefore, coffee drinking is also associated with a decrease in the probability of having heart disease. This association appears to be significant.
• Fast Food Spend: A one-unit increase in fast food spend is associated with a 0.0024 increase in the log odds of having heart disease. This association appears to be heavily significant.
• Income: Based on a p-value of 0.8182, we conclude that income is not significantly associated with heart disease.

Name: Biopsy Data on Breast Cancer Patients
Documentation: See reference on the data downloaded and conditioned for R here https://cran.r-project.org/web/packages/MASS/MASS.pdf
Sampling details: This breast cancer database was obtained from the University of Wisconsin Hospitals, Madison from Dr. William H. Wolberg. He assessed biopsies of breast tumours for 699 patients up to 15 July 1992; each of nine attributes has been scored on a scale of 1 to 10, and the outcome is also known. The dataset contains the original Wisconsin breast cancer data with 699 observations on 11 variables.

The biopsy data contains 699 observations of 9 continuous variables, V1, V2, …, V9.

The dataset also contains a character variable: id, and a factor variable: class, with two levels (“benign” and “malignant”).

# check variable types
str(biopsy)

'data.frame':   699 obs. of  11 variables:
 $ ID   : chr  "1000025" "1002945" "1015425" "1016277" ...
 $ V1   : int  5 5 3 6 4 8 1 2 2 4 ...
 $ V2   : int  1 4 1 8 1 10 1 1 1 2 ...
 $ V3   : int  1 4 1 8 1 10 1 2 1 1 ...
 $ V4   : int  1 5 1 1 3 8 1 1 1 1 ...
 $ V5   : int  2 7 2 3 2 7 2 2 2 2 ...
 $ V6   : int  1 10 2 4 1 10 10 1 1 1 ...
 $ V7   : int  3 3 3 3 3 9 3 3 1 2 ...
 $ V8   : int  1 2 1 7 1 7 1 1 1 1 ...
 $ V9   : int  1 1 1 1 1 1 1 1 5 1 ...
 $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...

There is one incomplete variable V6 = “bare nuclei” with 16 missing values.

• remember the package skimr for exploring a dataframe?

# check if vars have missing values
biopsy %>% 
  # select only variables starting with "V"
  skimr::skim(starts_with("V")) %>%
  dplyr::select(skim_variable, 
                n_missing)

# A tibble: 9 × 2
  skim_variable n_missing
  <chr>             <int>
1 V1                    0
2 V2                    0
3 V3                    0
4 V4                    0
5 V5                    0
6 V6                   16
7 V7                    0
8 V8                    0
9 V9                    0

We can decide what to do in these cases (informed by our knowledge of the dataset):

• Option 1) We drop the observation with incomplete data (i.e. with missing values for V6 = “bare nuclei”) with 16 missing values.

# remove rows with missing values
biopsy_drop <- biopsy %>% 
  dplyr::filter(!is.na(V6))

mean(biopsy_drop$V6)

[1] 3.544656

• Option 2) We impute the missing values with the mean of the variable V6 = “bare nuclei”.

# impute missing values with the median of the variable
biopsy_impute <- biopsy %>% 
  dplyr::mutate(V6 = ifelse(is.na(V6), median(V6, na.rm = TRUE), V6))

mean(biopsy_impute$V6)

[1] 3.486409

• Biopsied cells of 700 breast cancer tumors, used to determine if the tumors were benign or malignant.
• This determination was based on 9 characteristics of the cells, ranked from 1(benign) to 10(malignant):
  • 1) Clump Thickness – How the cells aggregate. If monolayered they are benign and if clumped on top of each other they are malignant
  • 2) Uniform Size – All cells of the same type should be the same size.
  • 3) Uniform Shape If cells vary in cell shape they could be malignant
  • 4) Marginal Adhesion – Healthy cells have a strong ability to stick together whereas cancerous cells do not
  • 5) Single Epithelial Size – If epithelial cells are not equal in size, it could be a sign of cancer
  • 6) Bare nuclei – If the nucleus of the cell is not surrounded by cytoplasm, the cell could be malignant
  • 7) Bland Chromatin – If the chromatin’s texture is coarse the cell could be malignant
  • 8) Normal Nucleoli – In a healthy cell the nucleoli is small and hard detect via imagery. Enlarged nucleoli could be a sign of cancer
  • 9) Mitosis – cells that multiply at an uncontrollable rate could be malignant

• The explanatory variable(s) (clump_thickness, …, mitosis) can be renamed for better readability
• The observations with missing values (bare_nuclei) are removed for simplicity
  
• Patient ID (id) can be dropped as it is not used in this analysis

# Create a clean version of the dataset
biopsy_clean <- biopsy %>%
  # rename the columns (new = old)
  rename(
    id                = ID,
    clump_thickness   =  V1,
    uniform_size      =  V2,
    uniform_shape     =  V3,
    marginal_adhesion =  V4,
    single_epith_size =  V5,
    bare_nuclei       =  V6,
    bland_chromatin   =  V7,
    normal_nuclei     =  V8,
    mitosis           =  V9,
    class             =  class) %>% 
  # remove rows with missing values
  na.omit(bare_nuclei) %>% 
  # remove the id column
  select(-id)

# check the structure of the dataset
paint::paint(biopsy_clean)

data.frame [683, 10]
clump_thickness   int 5 5 3 6 4 8
uniform_size      int 1 4 1 8 1 10
uniform_shape     int 1 4 1 8 1 10
marginal_adhesion int 1 5 1 1 3 8
single_epith_size int 2 7 2 3 2 7
bare_nuclei       int 1 10 2 4 1 10
bland_chromatin   int 3 3 3 3 3 9
normal_nuclei     int 1 2 1 7 1 7
mitosis           int 1 1 1 1 1 1
class             fct benign benign benign benign benign ma~

• We fit a logistic regression model to the biopsy_train data using:
  • the glm function with argument family = binomial to specify the logistic regression model;
  • and with Class ~ . to specify an initial model that uses all the variables as predictors (backward elimination approach).

# Building initial model 
model = stats::glm(class ~ . , family = binomial, data=biopsy_train)

• Table 2 shows the model summary, with the coefficient estimate for each predictor.
  • the broom::tidy function converts the model summary into a data frame.

broom::tidy(model) %>% 
  mutate('Sign.lev' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  knitr::kable() 

In logistic regression, the coefficients are in the form of the natural logarithm of the odds of the response event happening (i.e. \(Y_i = 1\)):

\[logit(p_i) = \ln\left( \frac{p_i}{1-p_i} \right) = -9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis\]

However, with some algebraic transformation, the logit function can be inverted to obtain the probability of the response event happening as a function of the predictors:

\[p_i = \frac{1}{1 + e^{-(-9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis)}}\]

This equation represents the logistic regression model’s best-fit line.

• The model includes all variables, but we could make it more parsimonious by removing variables that are not significant!
• We use a statistic called the Akaike Information Criterion (AIC) to compare models.
  • AIC’s calculation gives a penalty for including additional variables.
  • The model with the lowest AIC is considered the best.

# For example let's fit a model without the variable `uniform_size`
model2 = glm(class~ .-uniform_size, family = binomial, data=biopsy_train)

According to the AIC values, the model2 seems better (AIC is lower).

# Compare the AIC values of the 2 models
tibble(Model = c("model", "model2"), 
       AIC = c(AIC(model), AIC(model2) )) %>% 
  kable()

• The MASS package’s function stepAIC enables to perform a systematic model selection (by AIC):
  • The direction argument specifies the direction of the stepwise search.
  • The trace argument (if set to TRUE) prints out all the steps.
• The best_model has removed these variables:
  • uniform_shape
  • single_epith_size
  • normal_nuclei
• The best_model has the lowest AIC value (from 100 to 98.5), despite a higher Residual Deviance than the full model (from 80 to 80.6), albeit by a very slight amount.

# Select the best model based on AIC
best_model <- MASS::stepAIC(model, direction = "both", trace = FALSE)

# Compare the AIC values of full and best model
tibble(Model = c("model", "best_model"), 
       AIC = c(AIC(model), AIC(model2)),
       Deviance = c(deviance(model), deviance(best_model))) %>% kable()

A 3D scatterplot of observations shows the first 3 principal components’ scores.

• For this one, we need the scatterplot3d() function of the scatterplot3d package;
• The color argument assigned to the Label variable;
• To add a legend, we use the legend() function and specify its coordinates via the xyz.convert() function.

# 3D scatterplot ...
plot_3d <- with(PC_scores, 
                scatterplot3d::scatterplot3d(PC_scores$PC1, 
                                             PC_scores$PC2, 
                                             PC_scores$PC3, 
                                             color = as.numeric(Label), 
                                             pch = 19, 
                                             main ="Figure 3: 3D Scatter Plot", 
                                             xlab="PC1",
                                             ylab="PC2",
                                             zlab="PC3"))

# ... + legend
legend(plot_3d$xyz.convert(0.5, 0.7, 0.5), 
       pch = 19, 
       yjust=-0.6,
       xjust=-0.9,
       legend = levels(PC_scores$Label), 
       col = seq_along(levels(PC_scores$Label)))

Biplots have two key elements: scores (the 2 axes) and loadings (the vectors). As in the scores plot, each point represents an observation projected in the space of principal components where:

• Biopsies of the same class are located closer to each other, which indicates that they have similar scores referred to the 2 main principal components;
• The loading vectors show strength and direction of association of original variables with new PC variables.

As expected from PCA, the single PC1 accounts for variance in almost all original variables, while V9 has the major projection along PC2.